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Abstract. We present a class of symplectic integrators adapted for the integration of perturbed Hamil- 
tonian systems of the form H — A + eB. We give a constructive proof that for all integer p, there exists 
an integrator with positive steps with a remainder of order 0(t p £ + t 2 s 2 ), where r is the stepsize of 
the integrator. The analytical expressions of the leading terms of the remainders are given at all orders. 
In many cases, a corrector step can be performed such that the remainder becomes 0(t p e + r 4 e 2 ). The 
performances of these integrators are compared for the simple pendulum and the planetary 3-Body problem 
of Sun- Jupiter-Saturn. 
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1. Introduction 

Symplectic integrators, due to their good stability properties are now currently used for 
long time integrations of the Solar System, starting with the work of Wisdom and Holman 
(1991). Despite some improvement resulting from a good choice of initial conditions (Saha 
and Tremaine, 1992) or a corrector to the output of the numerical integration (Wisdom 
et al, 1996), it is surprising that the integration method which is currently used in most 
computations (see Duncan et al, 1998) is the celebrated 'leapfrog' method of order 2 
(Ruth, 1983). A reason for this choice is probably due to the fact that the methods of 
higher order which have been found by Forest and Ruth (1990) or Yoshida (1990) do not 
present very good stability properties for large stepsize, due to the presence of negative 
steps. 

In the present work, by considering perturbed Hamiltonians on the form H = A + 
eB were both A and B are integrable, we proove the existence of a class of symplectic 
integrators with positive steps which improve the precision of the integration by several 
order of magnitude with respect to the commonly used leapfrog method, and which present 
good stability properties at large stepsize. 

2. Lie formalism 

According to Yoshida, (1990), the search of symplectic integrators using Lie formalism 
was introduced by Neri (1988). Since, it was largely developped by Yoshida (1990), Suzuki 
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(1991, 1992), Koseleff (1993, 1996), and Mclachlan (1995, 1998). Let H(p,q) be an Hamil- 
tonian denned on R n x T n , where (p, q) are the actions and angle variables. Hamilton 
equations are 

dpj dH dqj dH 
dt dqj ' dt dpj 
and the Poisson bracket of /, g is defined on R n x T n by 

{/■«>- £££-2^ <*> 

j d Pj <% dqj d Pj 

If we denote x = (p, q), we obtain 

dx 

— = {H,x} = L H x. (3) 
dt 

where Ljj is the differential operator defined by L x f = /}. The solution x(t) of (3) 
with x(0) = xo is obtained formally as 

x(t) = T-L n H x = e tL »x . (4) 

A symplectic scheme for integrating (3) from t to t + r consists to approximate in a 
symplectic way the operator q tLh . Indeed, as H = A + eB, the Campbell-Baker-Hausdorf 
(CBH) theorem ensures that 

e rL H = e rL Ae rL eB + ^ _ ^ 

The operator S± = e TL A e rL sB pj-gyi^gg the most simple symplectic sheme for such 
Hamiltonians. This can be generalized with a combination of several steps involving suc- 
cessively A and eB in order to obtain integrators of higher orders. A general integrator 
with n steps will be 

where the constants (q, di) will be chosen in order to improve the order of the integrator. 
Using CBH theorem, and the linearity of the Lie derivative, we are ensured of the existence 
of a formal series 

K = ki^A + ek h2 B + rek 2 ,i{A, B} 

+T 2 ek 3A {A, {A B}} + T 2 e 2 h, 2 {{A, B}, B} 

+r 3 ek 4)1 {A, {A, {A, B}}} + r 3 s 2 k 4>2 {A, {{A, B}, B}} 

+r 3 e 3 k 4 , 3 {{{A, B}, B}, B} + 0(r 4 ) 

where the coefficients kij are polynomials in the (c m ,d n ), with rational coefficients, such 
that 

S n ( T ) = e ^L Ae dirL £B . . . e c n rL Ae d n rL £B = & tL k (g) 
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It should be noted that in order to define these expressions in a non ambigous way, one 
needs to decompose the Poisson brackets involving A and B over a basis of canonical 
elements of the free Lie algebra C(A, B) generated by A and B and the Poisson backet 
{ , }. Following Koseleff (1993), this is done here by using the Lyndon basis. The scheme 
S n {r) integrates in an exact manner the formal Hamiltonian K. A symplectic integrator 
for H = A + eB will be obtained at order p if K = A + sB + 0(t p ). In the most general 
way, this will be achieved by solving the algebraic equations 

fci,i = 1 ; fci )2 = 1 ; 
kij = for (i<p). (9) 

In particular, we have k± i = c\ + c 2 + • ■ ■ Cn = 1, k± 2 = d\ + <i 2 + ■ ■ ■ d n = 1, for p > 1. 



3. Symmetric integrators 

We will now restrict ourselves to symmetric integrators, that is integrators S n (r) such that 
Snir)- 1 = S n (-r). We will have 

-r Ljf( T ) = -r L^(_ T ) (10) 

thus if (— r) = K(t), and the formal Hamiltonian if (r) is even. As we distinguish A and 
e.B, we will have several classes SABAk and SBABk of symmetric symplectic operators 
defined by their prototypes 

SABA211 ' e CirLA e dlTLeB • • • Q dnTL ^B Q C n+lTL A ^d n TL eB _ _ _ tjl\TL eB e cirL A 

SABA2n+i '■ e ClTLA e dirIjEB • • • e Cn+1TLA e dn+1TLeB e Cn+irLA ■ ■ ■ e dlTLeB e ClTLA 
SBAB211 ' e dirLsB e C2TLA e d2TLeB • • • e rfnTi£S e Cn+1TiA e rfnTieS • • • e C2T ^ A e rflT ^ £S 

SBAB2n+l ' e dlTLsB Q C2tLa ■ ■ . g c n+ 1 ri A g 4+ 1 Ti E B g C„ + l tLa . . . QC2TL A ^dlTL eB 

The index of the integrator is the number of evaluations of A and B which are necessary 
for each step. With these notations, the classical leapfrog integrator can be considered as 
SBAB! = e ^ L - B e rLA ei L - B £ SBAB\ or as SAB Ai = e i LA e TL * B ei LA £ SABAi. In both 
cases, the integrator is of order 2 and the formal Hamiltonian is K = A + eB + 0(r 2 e). 
The fourth order solution found by Forest and Ruth (1990) or in an other way by Yoshida 
(1990) is either of the form SAB A3 or SBAB 3 that is, for SBAB 3 

SFRA3 = e dirL i B e c 2 TL A e d2rL EB e c 3 TL Ae d 2 TL eB ^itL a ^tL cB ^y£^ 

with 

'eg + 2c 2 = 1 

di + d 2 = 1/2 (13) 
1/12 - 1/2 c 2 + 1/2 c 2 2 + c 2 di - c 2 2 di = 0; V ; 

L -1/24 + 1/4 c 2 - c 2 di + c 2 d\ = 
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This system has a single real solution with approximate values d\ « 0.6756, C2 ~ 1.3512, 
d,2 ~ —0.1756, C3 ~ —1.7024. The problem with this integrator, is that due to the presence 
of negative time steps, the absolute value of the time steps remains high, and for large 
stepsizes, at an equivalent cost, the leapfrog integrators becomes more effective. In fact, 
Suzuki (1991) has demonstrated that it is not possible to obtain integrators of order 
p > 2 with only positive steps. The problem of the negative stepsize can nevertheless be 
overcome. 



4. Integrators for perturbed Hamiltonian 

In the previous sections, we have not yet taken into account the existence of the small 
parameter e. Indeed, the terms of second order of K (7) are r 2 ek 3 i{A, {A, B}} and 
T 2 e 2 k3 : 2{{A, B}, B} which are respectively of order t 2 e and t 2 e 2 . One can thus try to 
cancel uniquely the largest term, that is ks,i = 0. This can be done using 

SABA 2 : e ClTiA e dirL£S e C2TLA e dlTLeS e ClTLA (14) 



or 



SBAB 2 : e dlTLsB e C2TLA e d2TL " B e CirLA e dlTL£B . (15) 
With the type SABA2, one obtains d\ = |, C2 = 1 — 2ci and 



KsabAz =A + eB +r 2 e{^ - -a + ^c?){A, {A, B}} 

+T 2 e 2 (-±+ 1 -c 1 ){{A,B},B} + 0(T*s) 



(16) 



As we search for only positive stepsize, we find a unique solution for cancelling the term 
in er 2 , that is 

c 2 = ^; ^ = \{i-^=)- dl = \- (17) 

with these coefficients, we obtain Ksaba 2 = A + eB + 0(r A e + t 2 e 2 ). In a similar way, we 
obtain the solution for SBAB2 

2 11 

&2 = - ; di = - ; c 2 = - ; (18) 

and as previously Ksbab 2 = A + eB + 0(r 4 e + t 2 e 2 ). Quite surprisingly, this latest 
integrator which is in most cases much more precise than the leapfrog integrator {SBAB\ ) 
at the same cost (see section 8), does not seem to have been used so far. 
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5. Higher orders 

It becomes then tempting to iterate this process at higher order. We will not try to remove 
the term of order r 2 e 2 , which is not the most important for large stepsize when e is small. 
We will search for solutions S n of the form SAB An or SBAB n for which the associated 
Hamiltonian Ks n verifies 

K Sn = A + eB + 0(r 2n e + tV) (19) 

For this, we need to cancell at all order p < 2n the coefficient k Pt \ of the single term of 
order r v e in the Lyndon decomposition of Kg n 

r p ek pA {A,{A,{A,...{A,B}}}...} (20) 

We thus need to compute the part of K$ n which is of degree < 1 in B. We will use some 
results on calculus on free Lie algebra for which the reader should refer to (Bourbaki, 
1972). We will call C(U, V) the free Lie algebra generated by U and V, endowed with its 
canonical associative structure. We will also use the symbol = for the equality in C(U, V) 
modulo terms of degree > 2 in V . We have the two lemmas (Bourbaki, 1972) 

LEMMA 1. 

e u V e ~u = e ad(U) y ( 21 ) 



where the exponential of X is formally defined as exp(X) = X^n^o X n /n\, and where the 
adjoint operator ad is defined as ad(X).Y = [X,Y]. 



1 _ e -ad(UY 



LEMMA 2. 

° u+v ^ v + * u {^uT) v - (22) 

The next result is a generalisation of a classical expansion at degree 1 in V of the Campbell- 
Baker-Haussdorff formula. 

PROPOSITION 1. Let 7 e R. Then there exists W G £{U, V) such that 

e 7U e v e (i-y)U = e w (23) 

with 

ad{U)^ adi ' U) , x 

W = U + e a d lu ) _ l V ( 24 ) 



that is 



Bp (7) 

W = U + V -^^ad(U) p V (25) 
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and where B n {x) are the Bernoulli polynomial defined as 



t e tx +00 t n 

= Y^B n (x)- (26) 
e — 1 ^— i n! 

n=0 

Indeed, the existence of W € C(U, V) satisfying the above relation is given by the CBH 
theorem, on the other hand, we have 

e iu e v e (i-^u ^ e u + e u e (n,-i)u Ve (i-y)u (27) 

and from lemma 1, this is also equal to 

e U + e U e ^-l)ad(U) v ^ (28) 

As for V = 0, we have W = U, we can set W = U + W\ , where W\ is of degree 1 in V, 
and from lemma 2 

(\ _ e -«i(tO\ 

^^h^rh < 29 > 

thus 

which ends the proof. For 7 = 1, we recover the CBH results. This result is then easily 
generalized to the case of multiple products. 

PROPOSITION 2. Let a, . . . c n , d 1: . . . , d n £ R, such that £? =1 a = 1. T/ien there exists 
W G £(t7, V) such that 

e cif/ e diy e c 2 [/ e d 2 y _ _ _ e c„i/ e d„v = e w ^ 

with 

W^ + E^ Urn : y (32) 



fc=l e 1 



w = ^ + E ( E * ^r) ad ^ p v ( 33 ) 



73 

p=0 \fc=l F 

mi/i 7 fc = ci + . . . + c fc . 

Dems. This is staightforward as soon as we remark that 



e aU e rfiV e c 2 u e d 2 v a c n U^d n V 



. . . e CnU e dnV = X d k e^ u Ve {1 ~ lk)u + e^ (34) 



k=i 
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Remark : As Bq(x) = 1, if Ylk=i ^fc = 1> we have 

B P {lk) 



w = u+v+j2\ 



p=l \fe=l 



ad(£/) p V . 



(35) 



6. Computation of the coeffcients 

Propostion 2 , applied with U = tLa and F = teLb, gives directly the algebraic equations 
which could then be solved for obtaining integrators of arbitrary order for perturbed 
systems. The problem is thus reduced to the search for coeffcients jk,dk such that 



Y. d *9{lk,t) = l+o(t N ) 



k=i 



for N as high as possible with 



g(x,t) 



te 



X t 



1 



(36) 



(37) 



That is, with J2k=i c k = 1) we wm have to solve an algebraic system of equations of the 
form 



^2 dk B (-f k ) = ^4 = 1 

k=l k=l 
n 

J2 d kB p {j k ) = forO<p<iV. 



(38) 



k=i 



It should be noted that all the integrators SABA n and SBAB n can be written on the 
general form (31) by taking d n = or c\ = in (31). Moreover, if we search for symmetric 
integrators, all the relations in (38) will be automatically fullfilled for odd values of p. 
In this case, we just have to consider even values of p, for which we give the Bernoulli 
polynomials up to p = 10. 



Bq(x 
B 2 (x 

B A {x 
B 6 (x 
B 8 (x 
B w (x) 



= 1 

1 2 
= --X+X 

o 

1 



hx 2 - 2x 3 

30 

1 x 2 5 x 4 

42 ~ Y + ~2~ " 

1 2x 2 7x 4 



3x 5 + x 6 



(39) 



30 + 3 



+ 



14 x b 



5_ 

66 



3x 2 



-Ax 1 + x b 
15 x 8 



+ 5x 4 -7x 6 + ^--5x 9 + x w 
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For example, the first integrators SABA2 = e clU e dlV e C2U e dlV e ClU will be obtained by 
solving the set of equations 

2ci + c 2 =1 

2d 1 = 1 (40) 

diB 2 (ji) + d 1 B 2 (^ 2 ) = 
with 71 = ci, 72 = ci + C2, thus 72 = 1 — 71 . As g(l — x,t) = g(x, —t), we have for all p 

B p (l-x) = (-l) p B p (x) (41) 

and the previous system reduces to 

rfi = 1/2 c 2 = l-2ci B 2 (ci)=0 (42) 

and we recover the previous results. For SAB A3 = e clC/ e (ilV e C2C/ e (i2 ^e C2f; e dl ^e clt/ we have 

ci + c 2 = 1/2 

d 2 + 2di =1 / 43 n 

diS 2 (7i) + d 2 B 2 (j 2 ) + ^1^2(73) =0 
diB 4 (7i) + ^2^4(72) + ^1^4(73) = 

with 71 = ci, 72 = ci + C2 = 1/2, and 73 = ci + C2 + C2 = 1 — ci. We have thus 
£2(72) = -1/12, ^4(72) = 7/240, #2(73) = # 2 (ci), B 4(73) = B 4 (ci). We are thus left with 

c 2 = 1/2 - ci 

d 2 = l- 2di / 44 x 
diS 2 (ci) - (1 - 2di)/24 = V 7 

di5 4 (ci) + 7(1 - 2di)/480 = 

The resolution of this system is made easily and provide a single solution for which all the 
coefficients Ci,di are positive 

5-VT5 s/15 , 5,4 

ci = ^-; c 2 = — ; * = -; d 2 = - (45) 

This can be continued at all orders, but algebraic equations becomes more complicated 
as the order increases. The symplectic integrators up to order 10 are listed in Table I. 
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Table I. Coefficients of the integrators SABA n and SBAB n up to n = 10. 



SABA! 



ci 1/2 | di 1 

SABA 2 


ci 1/2 -a/3/6 
c 2 \/3/3 


cfi 1/2 




ci 1/2-^15/10 
c 2 \/l5/10 


di 5/18 
d 2 4/9 




ci 1/2 - \/525 + 70\/30 /70 

c 2 (\/525 + 70\/30 - \/525 - 70\/3o) /70 

c 3 \/525 - 70\/30 /35 


di l/4-\/30/72 
d 2 l/4 + x/30/72 




ci 1/2 - (\/490 + 42x/l05 + \/490 - 42x/T05)/84 

c 2 \/490 - 42\/T05/42 

c 3 (\/490 + 42x/T05 - \/490 - 42x/T05)/84 


di (322 - 13\/70)/1800 
d 2 (322 + 13\/70)/1800 
d 3 64/225 




ci 0.033765242898423986093849222753002695 
c 2 0.135630063868443757075450979737044631 
c 3 0.211295100191533802515448936669596706 
c 4 0.238619186083196908630501721680711935 


di 0.085662246189585172520148071086366447 
d 2 0. 180380786524069303784916756918858056 
d 3 0.233956967286345523694935171994775497 




ci 0.025446043828620737736905157976074369 
c 2 0. 103788363371682042331 162455383531428 
c 3 0.167843017110998636478629180601913472 
c 4 0.202922575688698583453303206038480732 


di 0.064742483084434846635305716339541009 
d 2 0.139852695744638333950733885711889791 
d 3 0.190915025252559472475184887744487567 
d 4 256/1225 



SABAs 



ci 0.019855071751231884158219565715263505 

c 2 0.081811689541954746046003466046821277 

c 3 0. 135567033748648876886907443643292044 

c 4 0.171048883710339590439131453414531184 

c 5 0.183434642495649804939476142360183981 



di 0.050614268145188129576265677154981095 

d 2 0.111190517226687235272177997213120442 

d 3 0.156853322938943643668981100993300657 

d 4 0.18134189168918099148257522463859781 



SABA g 



ci 0.015919880246186955082211898548163565 

c 2 0.066064566090495147768073207416968997 

c 3 0.1 1 1329837313022698495363874364130346 

c 4 0. 144559004648390734135082012349068788 

c 5 0.162126711701904464519269007321668304 



di 

rf2 

d 3 
d 5 



0.040637194180787205985946079055261825 
0.090324080347428702029236015621456405 
0.130305348201467731159371434709316425 
0.156173538520001420034315203292221833 
16384/99225 



ci 0.013046735741414139961017993957773973 

c 2 0.054421580914093604672933661830479502 

c 3 0.092826899194980052248884661654309736 

c 4 0.123007087084888607717530710974544707 

c 5 0.142260527573807989957219971018032089 

c 6 0.148874338981631210884826001129719985 



di 0.033335672154344068796784404946665896 

d 2 0.074725674575290296572888169828848666 

d 3 0.109543181257991021997767467114081596 

d 4 0.134633359654998177545613460784734677 

d 5 0.147762112357376435086946497325669165 
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Table I. 



SBABi 



c 2 1 | di 1/2 

SBAB 2 


C2 1/2 


di 1/6 
d 2 2/3 


SBABz 


c 2 1/2-^/5/10 
c 3 VE/5 


di 1/12 
d 2 5/12 


SBAB 4 


c 2 l/2-V^77/2 


di 1/20 
d 2 49/180 
d 3 16/45 




c 2 1/2-^3 + 6/^7/6 

cs (^/3 + 6/V7-y / 3-6/\/7)/6 

c 4 ^1/3-2/3^ 


di 1/30 

d 2 (14-\/7)/60 

d 3 (14 + \/7)/60 




c 2 1/2- y / (15 + 2yi5)/33/2 
c 3 ^5/22- \/5/33/2 
c 4 ^5/44-^573/22 


di 1/42 
d 2 31/175 - VW5/20 
d 3 31/175 + VW5/20 
d 4 128/525 




c 2 0.064129925745196692331277119389668281 
c 3 0.140019983538232156596467514911355124 
c 4 0.191200481765331716687926735526300967 
c 5 0.209299217902478868768657260345351255 


di 1/56 
d 2 0.105352113571753019691496032887878162 
d 3 0.170561346241752182382120338553874086 
d 4 0.206229397329351940783526485701 104895 




c 2 0.050121002294269921343827377790831021 
c 3 0.111285857950361201933229908663497754 
c 4 0.157034407842279797367566679191341619 
c 5 0. 181558731913089079355376034354329607 


di 1/72 
d 2 0.082747680780402762523169860014604153 
d 3 0. 137269356250080867640352809289686363 
d 4 0.173214255486523172557565766069859144 
d 5 2048/11025 


SBAB g 


c 2 0.040233045916770593085533669588830933 
c 3 0.090380021530476869412913242981253705 
c 4 0. 130424457647530289670965541064286364 
c 5 0. 156322996072028735517477663386542232 
c 6 0.165278957666387024626219765958173533 


di 1/90 
d 2 0.0666529954255350555631 13585377696449 
d 3 0.112444671031563226059728910865523921 
d 4 0.146021341839841878937791128687221946 
d 5 0. 163769880591948728328255263958446572 


SBAB W 


c 2 0.032999284795970432833862931950308183 
c 3 0.074758978372457357854928159995462766 
c 4 0. 109624073333469706075726923315353220 
c 5 0. 134738595704632807519526226959347078 
c 6 0. 147879067793469695715955757779528754 


di 1/110 
d 2 0.054806136633497432230701724790175355 
d 3 0.093584940890152602054070760949717460 
d 4 0.124024052132014157020042433210936377 
d 5 0. 14343956238950404433961 1201665767616 
d 6 32768/218295 
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7. McLachlan solution 



While we were writing a first version of this work, we realized that McLachlan (1995) 
had already found all the previous integrators. The paper of McLachlan is obviously 
not well-known to astronomers, otherwise they would have used at least the integrators 
SABA2, SBAB2, SAB A3 and SBAB% which have very good properties 1 . McLachlan just 
makes the computations on a very simple example for which the integration of the equa- 
tions reduces to a simple integral. He then claims that this is representative of the most 
general case. Although this may be true, the argument is not as straightforward as the 
constructive method which is presented here. On the other hand, the final remarks of 
McLachlan (1995) can be adapted here to complete the present proof and to provide the 
expression for the coefficients of these symplectic integrators at any order. Indeed, if we 
observe that 



j _ 



t Jo 



e xt dx , (46) 



and that (e* — l)/t = O(l), the problem of finding dk,jk verifying (36) is equivalent to the 
search of weights dk and nodes 7^ such that 

n ,.\ 

V4e 7fct = / e xt dx + o(t N ) (47) 

The solution of this problem is known classically as the Gauss integration formula. The 
values of 7/% are given by 7^ = (1 + x^)/2 where x k are the roots of the degree n Legendre 
polynomial P n {x). The the associated weights dk are all positive and are given by 

d k = o (48) 

{l-x 2 k ){P' n {x k )) 2 

More precisely, if we consider an integrator of type 

SABA n : e ClU e dlV e C2U e d2V . . . e CnU e d " v e c "+ lU , (49) 

without any assumption of symmetry, we will have, in the above formula d n+ i = 0, thus, 
for k = 1, . . . ,n, the coefficients 7^ = (l + Xk)/2 where x& are the roots of P n (x). All Xk are 
in the interval [—1, 1]. We will thus have 7^ G [0, 1]. If we put the 7^ in ascending order, 
the values of the coefficients Ck = 7fe+i — 7fe are all positive and c n+ \ = 1 — 7„. Moreover, 
the roots of the Legendre polynomial are symmetric with respect to zero. The 7^ are thus 
symmetric with respect to 1/2 and so will be the c& and dk- The symplectic integrator is 
thus symmetric, and this hypothesis was not necessary. This is not the case for Lie algebra 



1 The first integrators of the family (SABA2 , SBAB2 , SAB A3 and SBAB3) have been also recently 
reported by Chambers and Murison (2000). The integrator SBAB2 is mentioned in the book of E. Forest 
(1998). 
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symbolic computation, where the assumption that the integrator is symmetric decreases 
in a large amount the number of the variables. For an integrator of type 

SBAB n : e dlV e C2U e d2V . . . e CnU e dnV e c "+ lU e d "+ lV , (50) 

we need to set 71 = c\ = in formula (47), which means that in the integration formula, 
one node is fixed to an extremity of the interval [0, 1]. On the other hand, we have 7 ra +i = 
Y^,i=i c i = 1- The problem is thus to find nodes and weights for a Gauss formula with 
fixed nodes at the boundary of the interval of integration. The solution is given by the 
Gauss-Lobatto formulas (Abramovitz and Stegund, 1965). For k = 2, . . . , n, we have = 
(1 + X] c )/2 where Xk are the n — 1 roots of Pn(x), and 

di = d n+1 = 1 ; d k = — — ) 2 ; for k = 2, . . . , n . (51) 

n(n + l) n{n + \){P n {x k )Y 

As previously, the integrators are symmetrical. These relations thus allow us to obtain 
in a straightforward manner symplectic integrators for perturbed systems at any order 
without the need to solve algebraic equations which are difficult to handle at large orders. 
Moreover, it provides a demonstration that this solutions exists at any order, with positive 
coefficients , c4 . 



8. Numerical examples 

In this section, we will test the efficiency of the family of integrators SABA n and SBAB n 
on a simple pendulum example and on a planetary problem. For the simple pendulum 

P 2 

# = y+£COS<? (52) 

we apply directly the previous computations with A = p 2 /2 and eB = ecosq. For each 
value of the stepsize r, we have measured the maximum difference between the energy at 
the origine and the computed energy along the trajectories, over a time T = 25000. This 
comparison is performed for e = 0.1 and e = 0.001 (Fig. 1). For SABAn or SBAB n , the 
logarithm of differences are plotted versus log(r'), where r' = r/n. In such a way, as n 
is the number of evaluations of exp^rL^) and exp(drLfi) for the given integrator, the 
integrators are compared at constant cost. As expected, for sufficiently small stepsize, the 
residuals behave as t 2 e 2 for n > 2, and as T 2 e for the leapfrog integrator (n = 1). It is 
also clear that for small stepsize, nothing is really gained by increasing the order of the 
integrator (n), beyond n = 2. 

For large stepsize, this is not true, as the term T n+2 e or r ra+3 e (see next section) is 
still dominant, and we observe an increase of the slope with the order of the integrator, 
until unstabilities appear, probably due to the divergence of the remainders (it should be 
reminded that if a stepsize of 1 is used for SABAi, a stepsize of n is used for SABA n ). 
In most cases, n = 3 or n = 4 seems to be the best choices. 
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Figure 1. Fig. 1-6. Logarithm of relative energy error plotted versus log(r'), where r' = r/n, r the stepsize, 
and n is the index of the method (and the curve) for the various symplectic integrators of the family SABA„ 
and SBABn- Fig. 1. Simple pendulum with e = 0.1 (a-b) and e = 0.001 (c-d). 



In the case of the planetary iV-Body problem, the situation is more complicated. The 
Hamiltonian is splitted in an integrable Keplerian part, A, and a perturbation, B, cor- 
responding to the mutual gravitational interaction of the planets. The Keplerian part is 
integrated in elliptical coordinates, while the perturbation (which is essentially a sum 
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Figure 2. Relative energy error versus stepsize for the Sun- Jupiter-Saturn problem in Jacobi coordinates 
for the family of integrators SAB An and SBAB n - 



of invert of the mutual distance of the planets) is integrated in rectangular cartesian 
coordinates (Wisdom and Holman, 1991). 

There are several possible choice of coordinates for this decomposition. The initial choice 
of Wisdom and Holman, (1991), was to use Jacobi coordinates. In this case, B is integrable, 
as it depends only on the positions q. In Poincare heliocentric coordinates (see Laskar and 
Robutel, 1995), the expressions are more simple, but the perturbation B needs to be 
splitted in two terms B = B\(p) + B2(q) which depends uniquely on the momentum p, or 
on the positions q. As the methods which are presented here depends only on the linear part 
(in Lb) of the integrator, they can be adapted in a straightforward manner to this case, by 
substituting in their expresions exp L b 1 exp L b 2 or exp L b 2 exp L b 1 to exp Lb ■ In doing so 
one needs to be sure that the final symplectic scheme is still symmetric, which will ensure 
that no terms of order 2 will appear in the decomposition of the corresponding formal 
Hamiltonian K in equation (7). The use of these coordinates for symplectic integrators 
was first proposed by Koseleff (1993, 1996). 

In the present case, we will use Jacobi coordinates, as this choice will be motivated by 
the next sections which require that B depends only on q. In Jacobi coordinates, we did the 
computation for the Sun- Jupiter-Saturn system over 25000 years (Fig. 2), and obtained 
very similar results as for the simple pendulum with e = 0.001. This is understandable as 
this is of the order of the ratio of perturbation due to the mutual interaction of the planets 
over the potential of the Sun. It can be clearly seen that for all n > 2, these integrators 
outperformed by several order of magnitude the precision of the leapfrog integrator, except 
for very large stepsizes. The best choices being again n = 3 or n = 4. In all figures, it is 
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very obvious that the r 2 e 2 term is the main limiting factor. We will now make an explicit 
computation of this term and present a strategy to get rid of it. 



9. Computation of the remainders 



We compute here the remainders of the symplectic integrators SABA n and SBAB n . By 
switching the role of U and V in (31), we obtain easily 

PROPOSITION 3. Let c 1 ,...c n ,d 1 ,...,d n £R, such that £? =1 a = £™=i <k = 1. Then 
there exists W 6 £(U, V) such that 

^u^v^u^v e c n u e d n v = e w (53) 

with 



-co / n 



W =U + V + [ d k J ad(Uy V 

P=Ak=1 P ' ' (54) 

with 5q = 0, 5k = d\ + . . . + dk, and where = is the equivalence modulo terms of degree > 2 
in U and V in C(U, V). 

If we apply this result to compute the largest term in the remainder of the previous 
symplectic integrators, we obtain for each integrator 



W = A + B + (± c k il^zA) {{A , B}, B}r^ 



\k=i 



\k=l P ' J 



(55) 



We can be more specific for the two classes of integrators SABAn and SBAB n by taking 
into account the fact that these integrators are reversible. In this case, each integrator of 
the classes SABA2 n , SABA2 n +i, SBAB 2n , SBAB 2n +i, with J27=i c i = Y^=i di = 1, is 
the time-r evolution of the flow of the Hamiltonian W, with the following remainders : 
- SABAin- w e have 2n + 1 steps with d 2 n+i = and, for p = 0, . . . , n 



C-n+l+p — Cn+l—p ] 1n+p — 1 Tn+1— p j 

dn+p = d n -\-i—p ; 5 n j r p =1 5 n —p ; 



(56) 
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which gives, after reduction of the symmetries, and d n = 1/2 — <5 n _i, c n+ \ = 1 — 2^ ri 
W = A + £B+ (S^Lb 2 (1/2) + CkB 2 (5 k -i)) {{A, B}, B}r 2 s 



2 

j2n+2 



2.2 

(57) 

+ ( 2 g dkBM \ J^yBr^e + G>(tV + r 2 ^e) 

- SABA2n+\'- We have 2n + 2 steps, with d 2n +2 = 0, 4+1 = 1 - 2<5 n , c„ + i = 1/2 - j n , 
and, for p = 0, . . . , n 



(58) 



Cn+l+p — c n+2— p j Tn+l+p — 1 7n+l— p i 

4+1+p 4+1— p j ^n+p 1 ^rt+l— p j 

which gives, after reduction of the symmetries 

W = A + eB + ^ c fc 5 2 (<5fc-i)j {{A 5}, £}r 2 e 2 

/ * « \ L 2n+4 (59) 

+ ( 4+1 B 2 n+4(l/2) + 2 £ 4 B 2n+4 (7fe) j ( 2w 4 +4)! ^ r2n+4£ 



+ 0(rV + r 2n+6 e) 



- SBAB 2 n'- This case is easily obtained by setting ci = in 5„4£M.2n+i- We obtain for the 
new Hamiltonian 



/n+1 



W = A + eB + [J2 c k B 2 (S k -i) {{A, B}, B}r 2 e 2 



\k=2 



n \ L 2n+2 (60) 

+ f 4+1 S 2n+2 (l/2) +2^4^+2(7^)1 {2n A +2)[ BT2n+2£ 

+0{t a s 2 + T 2n+A e) 
SBAB 2n+ i: This case is easily obtained by setting c\ = in SABA 2n + 2 . 



/ n+l \ 

W = A + eB + ( ^B 2 (l/2) + £ CjfcJ B 2 (4-i) J {{A, B}, B}rV 

(Ji+1 \ r 2n+4 

2 ^ 4^+4(7.) j -^yy Br 2 ^s + 0(r'e 2 + r 2 " +6 e) 



(61) 
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10. Correctors 

The integrators SABA n and SBAB n have very good properties for small values of the 
parameter e. Their numerical properties were studied in section 8. We have seen that the 
main limiting factor is the term {{A,B},B}, which order is r 2 e 2 . It would be of course 
very nice to get rid also of this term, but the result of Suzuki (1991) tells us that it is 
not possible to get rid simultaneously of the two terms {{A, B}, B} and {^4, {A, B}} with 
integrators having only positive values for the q, di constants. It is not forbidden to have 
negative values for some of the constants, but as J2 c i = J2^i = 1, having only positive 
constants ensures that the values of the constants becomes smaller as the order of the 
integrator increases. This prevents explosion of the coefficients of the remainders which 
are polynomial in the Ci,di. 

In order to get rid of the {{A, B},B} term, one can use an alternate strategy, which 
is possible when A is quadratic in the actions p, and B depends only on the positions q 
(this is in particular the case for the pendulum Hamiltonian, or for the iV-Body problem 
when expressed in Jacobi coordinates). In this case, {{A,B},B} depends only on q and 
is thus integrable. It is then possible to compute it, and to add an additional step to the 
integrator S of the form 

S c = e- T / 2cL {{A,B},B}g e -r/2cL {{A:B}tB} 

where c is the coefficient of {{A,B},B} in W (Eq. 57-61). The new corrected integrator 
Sc is still symmetric, and thus additional terms will appear only at order r 4 . The values of 
the coefficients c used in the correctors up to order 10 are listed in Table II. For some of the 
lowest orders, algebraic formulas can be given, but they become very rapidly cumbersome, 
and a better accuracy will be obtained by using the decimal value which is given here with 
40 digits. 



Table II. Coefficients for the correctors up to order 10 for SABAn and SBAB-, 



n 


CSABA„ 


CSBAB„ 


1 


1/12 


-1/24 


2 


(2 - V3)/24 


1/72 


3 


(54- 13V5)/648 


(13- 5V5)/288 


4 


0.003396775048208601331532157783492144 


(3861 - 791 v / 2T)/64800 


5 


0.002270543121419264819434955050039130 


0.002381486672953634187470386232181453 


6 


0.001624459841624282521452258512463608 


0.001681346512091906326563693215296434 


7 


0.001219643912760418472579211822331645 


0.001251765616039400003072516100251191 


8 


0.000949308177745602234792177503535054 


0.000968797968073688571654684208462982 


9 


0.000759846022860436646358196674176815 


0.000772349023999952078227686810260323 


10 


0.000621934331486166426497049845358646 


0.000630320044163167840798638762665112 
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The plots of the residuals for these new integrators are presented in the case of the 
pendulum with e = 0.1 and e = 0.001 (Fig. 3-4), and the Sun- Jupiter-Saturn problem 
in Jacobi coordinates (Fig. 5). As we attain now the limitation due to round-off errors, 
computations were performed also in quadruple precision. It is clear that now the slope of 
the residuals corresponds to the r 4 terms and that we got rid of the r 2 e 2 term. 
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Figure 3. Relative energy error versus stepsize for the simple pendulum with s = 0.1 for SABA n and 
SBABn with correctors. 



11. Composition of integrators 

The corrector method of section 10 provide a family of integrators SABAcn, SBABcn of 
order 4 in r and higher order in e with remainders 0(r 4 e 2 ) + 0{r k e) with k = n + 2 for n 
even, and k = n + 3 for n odd. These integrators have very good numerical properties, but 
it is still possible to improve them by using the composition method of Yoshida (1990). 
Indeed, if Sir) is an integrator of order 2k, then it is possible to find c such that 

S(t)S(ct)S(t) (63) 

is an integrator of order 2k + 2. Indeed, the symmetry of the integrator ensures that there 
are no terms in r 2k+1 in the remainders, and a straightforward computation gives the 
condition of cancellation of the terms in r 2fc 

c 2fc+i + 2 = o (64) 

i 

that is c = — 2 2k + 1 . It should be noted that as c is close to —1, the cost of this composition 
scheme, which we will denote S 2 , is roughly 3 times more expensive than the initial 
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Figure 4- Relative energy error versus stepsize for the simple pendulum with e = 0.001 for SABAn and 
SBABn with correctors in double (a-b) and quadruple (c-d) precision. 

integrator S. Practically, we do one step forward, one step backward, and then one step 
forward again. Nevertheless, if one generalises this sheme to a composition S 2m defined as 



5 2m (r) = S m {r)S(cT)S m {r) 



(65) 



the condition (64) gives c = —(2m) 2fc + 1 . Usually c is still not very large, and the additionnal 
backward step becomes negligeable for large values of m. Unfortunately, as one would 



20 



J. Laskar and P. Robutcl 




expect, when m increases, the size of the remainders also increases and when we analyse 
the precision versus cost, it appears that we gain only for small values of m (Fig. 6). These 
integrators are still interesting, especially when one searches for high accuracy, which 
means small step size. 
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Figure 6. Relative energy error versus stepsize for the simple pendulum with e = 0.001 for the composition 
of SAB An and SBABn for n = 2 (a), n = 3 (b), n = 4 (c), and n = 5 (d). The index of the curve 
corresponds to the number of iterates 2m in the composition method (Eq.65). 
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12. Miscellanous remarks 

12.1. Integrals 

The following result is obtained immediately: 

PROPOSITION 4. Let H = A + B. If F is an integral of H and F commutes with A 
({A, F} = 0), then F is a true integral of the symplectic integration of H by any any of 
the integrators constructed above. 

Indeed, as {A, F} = 0, and {H,F} = 0, we have {B,F} = 0, and thus F commutes 
with any element of the free Lie algebra £(A, B). Thus, if the integrator <S(r) is defined 
by <S(t) = e rLw where W € £(A,B), we have {W, F} = 0. In particular, in Jacobi 
or heliocentric coordinates, the angular momentum depends only on the action variables 
and thus commute with the Hamiltonian of the Keplerian unperturbed problem. The 
angular momentum is thus an exact integral of the symplectic integration of the iV-body 
problem. In constrast, the initial Hamiltonian is only an approximate integral (at order 
0(T p e 2 ) + 0(T k e)). This feature can be used to check for the accumulation of errors in the 
integration. 

12.2. Non Hamiltonian systems 

In fact, the present results apply to general first order differential equations, and not only 
for Hamiltonian systems. Indeed, the only properties which are used are formal properties 
of the Lie algebra of the Lie derivatives along the vector fields defined by A and B. If a 
differential system of order 1 can be written on the form 

X = (L A + L B )X (66) 

where La and Lb are differential operators, for which the two systems X = LaX and 
X = L B X are integrable, then the symplectic integrators defined above will apply in the 
same way. Even more, if F is an integral of the system (66) such that La F = and 
Lb F = 0, then F is also an integral for the symplectic integrator. 

12.3. H = A + B 1 + B 2 

It happens very often that the perturbation is not integrable, but can be splitted in two 
parts B = B1 + B2 which are integrable separately (this is the case in Poincare heliocentric 
coordinates). As was already stated, the integrators SAB A n and SBAB n can be used 
provided some small modifications, but it will not be possible to use the correctors as 
defined in section 10. 
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We have presented here a new and constructive proof for the existence at all orders of 
the families of symplectic integrators SABA n and SBAB n , which were first described 
by McLachlan (1995). We have also obtained the expressions of the leading terms of the 
remainders for all n. These integrators are particularly adapted to perturbed Hamiltonian 
systems of the form H = A+eB, where A and B are integrable separately, and in particular 
for planetary iV-body problems. 

Moreover, when A is quadratic in the actions p and B depends only on the positions q, 
the new family of integrators SABAcn an d SBABcn given in section 10 provide integration 
scheme which is of order 4 in r, and has a remainder of the order of 0(r 4 e 2 + T p e), where 
p = n + 2orp = n + 3. For practical use, it seems that the integrators for n = 3 or n = 4 
are the most efficients. Despite they require additional computations for the corrector, 
the corrected integrators SABAcn an d SBABcn will beat the simple integrators SABAn 
and SBAB n in many occasions, but unless one search for very high accuracy with small 
stepsize, composition as described in section 11 is usually not very useful. 

All the integrators which are presented here have only positive stepsize, except for the 
corrector. It should still be investigated whether some integrators of order 4 with negative 
stepsize could be useful. 
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